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ABSTRACT 

Magnetorotational turbulence and magnetically driven disc winds are often considered as separate processes. However, realistic as- 
trophysical discs are expected to be subject to both effects, although possibly at different times and locations. We investigate here the 
potential link between these two phenomena using a mixed numerical and analytical approach. We show in particular that large-scale 
MRI modes which dominate strongly magnetised discs (plasma ji ~ 10) naturally produce magnetically driven outflows in the non- 
linear regime. We show that these outflows share many similarities with local and global disc wind solutions found in the literature. 
We also investigate the 3D stability of these outflows and show that they are unstable on dynamical timescales. The implications of 
these results for the transition between a jet-emitting disc and a standard "viscous" disc are discussed. 

Key words, magnetohydrodynamics (MHD) - instabilities - jets and outflows 



1. Introduction 

Accretion discs are found around several kind of astrophysical 
objects: from young stars (protoplanetary discs) to supermassive 
black holes in AGNs. The dynamics of these discs is however 
poorly understood. It is known that, on average, matter moves 
inward resulting in the accretion of mate rial onto the central ob- 
ject. Dynamically, one can show easily dLvnden-Bell & P ringle 
1974) that this accretion is controlled by the loss of angular mo- 
mentum of the falling gas. Hence, if one is to predict the accre- 
tion rate and survival time of discs, one needs to understand the 
angular momentum dynamics of these objects. 

For accretion to happen on timescales compatible with ob- 
servations, angular momentum must be tr ansported rather effi- 
cientl y. This has led to the a disc model (Sh akura & Sunvaevl 
1 1973b in which angular momentum is transported radially out- 
ward in the disc by a prescribed "turbulent viscosity". The 
origin of this turbulent viscosity has been debated for several 
decades and remains eve n today rather unclear. The magnetoro- 
tation al instability (MRI: Balbus & Hawlevlll99ll: lHawlev et~aTI 
1 1 1995b appears as the most promising way of producing turbu- 
lence in discs, though other processe s might also be at work 
in some regions (see lArmitage1l201 ll for a review). However, 
transporting angular momentum in discs is not the only way to 
cause accretion. Instead, one can suppose that angular momen- 
tum is "extracted" from a disc by a l arge-scale magnetic struc- 
ture anchored in the dis c midplane (jBlandford & Payne 1982; 
Pudri tz & Norman1 [l983). This sort of mechanism usually pro- 
duces magnetically driven outflows, also known as disc winds. 

It is most likely that these two mechanisms (angular momen- 
tum transport through turbulence and extraction by winds) actu- 
ally coexist in astrophysical objects. However, this link is, from 
a theoretical point of view, poorly understood. The first reason 
is that stati onary disc wind models re quire a strong magnetisa- 
tion p ~ 1 dFerreira & Pelletienl 19951) which quenches the MRI 



by magnetic tension effects (Bal bus & Hawlevlll99lh . At first 
sight, the conditions of existence for MRI turbulence and disc 
winds are therefore mutually exclusive. Moreover, a numerical 
computation including both a large-scale wind and small-scale 
turbulence in the disc midplane is technically very challenging: 
resolving turbulence in the bulk of the disc requires a large res- 
olution in the disc midplane and fully 3D simulations, whereas 
computing the wind itself requires a very large computational 
domain. Several authors have therefore computed disk winds us- 
ing a prescription to take into account the effects of disc turbu- 
lence (ess entially a turbulent visc osity and resistivity ), both an- 
alytic ally (Cas se & Ferre ira 2000) and numerically (Zann i et"a"D 
120071) . On the MRI side, most of our know l edge comes from 
shearing box simulations (lHawlev et aD 1995b IStone et alii 19961 : 
lLongaretti & LesuJ 1201 Oh . althou gh several au t hors have also 
considered global config urations (Hawley 2000; Froman g et al.1 
201 1; iFlock et alj|201 lb . These global configurations are how- 
ever limited to situations without any poloidal magnetic structure 
(or very weak mean poloidal magnetic fields) which precludes 
the production of outflows. 

More recently the production of outflows from MRI turbu- 
lence in shearing boxes has been studied in the limit of weak 
poloidal fields, /3 = 10 5 - 10 4 dSuzuki & Inutsukall2009l). A 
magn etocentrifugal mechanism similar to Blandford & Payne 
dl982l) has been identified in these simulations, with non-steady 
outflows starting abo ut two scale heights above the midplane 
dFromang et al.ll20l2l) . However, the outflow mass loss rate was 
shown to depend strongly on the box vertical size and aspect ra- 
tio. Moreover, the outflow for such a weak mean field is dynam- 
ically insignificant for the disc as it extracts a very small amount 
of angular momentum. In the stron g poloidal field regime Q3 < 
1) which is stable to MRI modes. lOgilviel d2012l) studied the 
production of quasi ID steady outflows in shearing boxes. This 
study demonstrated that some properties of global solutions were 
recovered in shearing boxes, although the procedure used in this 
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work was not used to look in the parameter regime unstable to 
the MRI. 

The aim of this paper is to study the potential link between 
the MRI and quasi-steady outflows which are found in global 
models and simulations. To this end, we study stratified shear- 
ing boxes threaded by a strong pol oidal field (B ~ 1 0) which are 
not far from the steady solutions of Qg ilviel (1201 2l) but are lying 
in the MRI unstable regime. We first present the shearing box 
model and the numerical methods we have used in this investiga- 
tion. We then look at the saturation of ID MRI modes which nat- 
urally produce outflows. These "MRI-outfiows" are compared 
to other solutions found in the literature, both local and global. 
We then study the 3D stability of these outflows and demon- 
strate that they are unstable on dynamical timescales. The con- 
sequences of this instability and its nature are briefly discussed. 
Finally, we summarise our findings and discuss their potential 
implications for astrophysical outflows and jets. 



2. Local model 

2.1. Equations 

MRI-related turbulence and the shearing box approximation 
have been extensively described in the literature. However, since 
the physics we are looking for involves outflows and mass losses, 
we recall here briefly the basic equations for the shearing bo x 
model. I nteres ted r eaders may also cons ult Hawlev et al.l (Il995l) . 
iBalbusI (l2003h and iRegev & Umurh an (2008) for an extensive 
discussion of this approximation. 

The shearing-box equations are found by considering a carte- 
sian box centred at a fiducial radius Ro, rotating with the disc at 
constant angular velocity Q = Q^(^o) and having dimensions 
(L x ,L y ,L z ) with Lj <s R$. We define r - Ro — > x, Rq4> — * y 
similarly to lHawlev et al.l (Il995l) . Furthermore, we will assume 
the disc follows an isothermal equation of state with a constant 
sound speed c. In this approximation, the MHD equations read: 



netic flux is conserved, we introduce the dimensionless magne- 
tization 



d,p + V -pu = 0, (1) 

d,pu + V • {pu ® u) = -c 2 Vp + (V X B) X B (2) 
-2pfl x u - pViff, 

d,B = V x (u x B), (3) 

where we have chosen the units so that po = 1 and iff is a local 
expansion of the effective gravitational potential around /?(>: 
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2.2 



iff = -qQ. l x A + -Q z z 



(4) 



in which we have considered a Keplerian disc having a rotation 
profile Qtf(V) oc f q with q = 3/2. 

One should note that the above set of equations admits a so- 
lution as a linear shear flow which is a local approximation to the 
Keplerian profile, u = -qQ.xe y . In the following, we will con- 
sider perturbations v (not necessarily small) to this Keplerian 
profile so that u — v - qQ.xe y . 

In the following, we consider a shearing box of size 
(L x ,L y ) = (16, 16), the unit of length being defined by the disc 
scale height: H — c/Q.. The time unit is Q. 1 and the velocity unit 
is c. We will assume p — 1 in the disc midplane at the hydrostatic 
equilibrium which sets our unit of mass. As shown below, we 
only consider the upper half of the disc, so that z e [0, zb] where 
Zb is the altitude of the upper z boundary condition. Unless oth- 
erwise stated, we assume zb = 6. Because the total vertical mag- 



SQc' 



(5) 



as a control parametetQ of our simulations, where 2 = 
J dV p/(L x Ly) is the equivalent surface density of the box. 

2.2. Numerical model 

2.2.1. Numerical method 

In order to investi gate the above syste m of equations, we use 
the PLUTO code dMignone et alJ feOO?). PLUTO is a finite vol- 
ume method using a Godunov scheme to integrate the equations 
in their conservative form. MH D terms are computed us ing the 
constrained transport method of lEvans & Hawlevl dl988l) which 
enforces V • B = at machine precision during the evolution 
of the physical system. We use the Roe method to solve the 
Riemann problem at cell boundaries. This choice was dictated 
by the presence of strong numerical instabilities with the HLLD 
solver when the plasma beta B = 2pc 2 /B 2 was too small (typi- 
cally B < 1). Moreover, in order to stabilise the code, we switch 
to an HLL solver (more diffusive) and a minmod slope limiter 
whenever the magnetisation is very large (typically B < 10" 4 ) 
in 3D runs. Only a few cells have such a small B and the pre- 
cise value of the threshold does not significantly change the out- 
come of the simulations. This however prevents the code from 
failing when the Alfven speed becomes too large. Similar tech- 
niques have been used in numerical studies of supersonic inter- 
stellar turbulence (Lem aster & Sto ne 2009). All the simulations 
discussed in this paper are summarised in Tab. [2] 

2.2.2. Boundary conditions 

In order to reduce the computational costs of the simulations, 
we compute only the upper half of the disc. The lower half 
is deduced by symmetry: p(-z) = p(z) ; v H (-z) = v H (z) ; 
v z (z) = -v z (z) ; B H (-z) = -B H (z) and B z (-z) = B z (z) where H 
stands for the horizontal (x, y) component of a vector. It should 
be noted that this symmetry is a natural symmetry of the equa- 
tions of motion. This implies that if the initial conditions satisfy 
this symmetry, the resulting solution will verify this symmetry 
at all times. 

The boundary conditions we impose are shear-periodic in the 
radial direction and periodic in the azimuthal direction. The mid- 
plane symmetry described above is imposed at z = 0. The upper 
boundary condition (z = Zb) is the most delicate part of the setup. 
Unless otherwise stated, we consider modified outflow boundary 
conditions where we enforce a strictly vertical poloidal field at 
the boundary: 



djKza) = d z v H (z B ) = d z B y (z B ) = 
B x (zb) = 



(6) 
(7) 



Surprisingly, strict outflow boundary conditions (zero gradient 
for all fields) prevent MHD driven winds to be produced. This 
point is discussed more extensively in 33.41 Boundary conditions 
used for each run are specified in Tab. [2] 

Because many of the simulations described here show an 
MHD-driven wind, a significant amount of mass is lost in our 



1 Note that this parameter is constant because we assume the total 
mass in the box to be constant 
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model. In order to mimic the mass inflow due to the material 
which would be accreted in a realistic disc including curvature 
effects, we have chosen to add mass in the midplane to com- 
pensate for wind losses. This is accomplished by adding mass at 
z < Zinj at each numerical time step. This procedure is in a sense 
a way to artificially introduce curvature effects in the shearing 
box, but it is clearly unsatisfactory. In particular, conservation of 
energy and momentum in the injection region z < Z; n j are not 
satisfied, which is clearly an issue if one considers the disc wind 
problem globally. In all the simulations discussed below all have 
Zinj =0.1, unless otherwise stated. Tests with Zi n j = 0.05 have 
shown that none of the result we discuss thereafter are signifi- 
cantly affected. 



2.2.3. Modified potential 

The vertical hydrostatic equilibrium described by <[2j> leads to a 
gaussian vertical density profile: 



p{z) = exp 



z 

Ztf 2 



Assuming a constant B, in the box we get an Alfven speed 
Va = B z l yfp oc exp(z 2 /4// 2 ) which increases very steeply as z 
increases. Because of the CFL condition, this leads to very small 
timesteps which dramatically increase the computational time. 
To reduce the computational costs, several of our simulations 
were performed using the modified potential: 



ip' = -qQ?x 2 + 



?Q 2 



1 - exp( 



-z 2 lzl) 



(9) 



This modified potential is roughly identical to the real potential 
for z < Zq but is less steep for z > Zo leading to a smoother den- 
sity profile in the hydrostatic equilibrium and therefore smaller 
Alfven velocities. It should be noticed that this problem arises 
only in the hydrostatic equilibrium without outflow. As we will 
show below, when an outflow is produced, the density profile is 
much smoother and we do not need a modified potential any- 
more. 

We have used the above modified potential with zo — 4 in 
our ID simulations to initiate the MHD wind in the tall box sim- 
ulations (z.b ^ 8). Once the quasi steady wind was formed, we 
relaxed back the potential to the original shearing box potential. 
Comparisons between the solutions obtained with and without 
the modified potential in the case zb — 6 have shown no differ- 
ence once the potential has been relaxed (runs IDRef — 1Dz6). 

3. One-dimensional MRI outflows 

In this section, we look at one-dimensional solutions of the equa- 
tions (Q]) — ([3]) i.e. (p(z, t), v(z, t), B(z, t)). This simplification al- 
lows us to isolate the physical mechanisms responsible for the 
MHD-driven outflows. 

We initialise our box with a disc in vertical hydrostatic equi- 
librium ([8]). We add a mean vertical field B : so that yt - 8 x 10~ 2 . 
In order to initialise the growth of MRI modes, we add a small 
perturbation B x - 0.02 sin(z) to the system. We show the tempo- 
ral evolution of this run (IDRef) in Fig.Q] 

3.1. From MRI modes to steady outflows 

At it can be seen from Fig. [T] we first observe the development 
of a linear MRI mode in the simulation b ox (t = 4.0). These lin- 
ear modes were described extensively bv iLatter et al.l (l2010h for 




Fig. 2. Growth rates of the largest stratified MRI eigenmodes as 
a function of the magnetizatio n parameter u. The se growth rates 
W were deduced from eq. (18) in lLatter et alJ d2010t) . 



stratified shearing boxes. The mode we observe in our particular 
setup is the n = 2 mode. This can be easily checked looking at 
th e shape of the per turbation and comparing to the eigenmodes 
of lLatter et all (l2010l) . Moreover, for fi = 8 x 10~ 2 , the n = 2 
mode is the most unstable mode in the system (Fig. [2j. Because 
the magnetic fluctuations are localised at z ~ H, the magnetic 
pressure tends to increase at that location, pushing toward the 
midplane the bulk of the disc and pushing up the disc atmo- 
sphere. At t — 7.0 - 8.0, the magnetic pressure is sufficiently 
large to push the atmosphere, creating a "bubble" of material. At 
this stage, the flow is no longer in a linear regime: the eigenmode 
is modified by nonlinear (e.g. magnetic pressure) terms. In the 
end, the system relaxes toward a quasi stationary state (f = 60). 
This implies that energy which is in injected by the MRI into the 
eigenmode is evacuated at the same rate in the outflow. 

It should be pointed out that we have totally ignored sec- 
ondary instabilities in this description. These secondary insta- 
bilities are usually x, y dependent modes which grow on the top 
of MRI eigenmode s (Good man & Xulll994t lLatter et al J 120091: 
iPessah &"G oodman 2009; Pessah l2010l) . It is widely believed 
that parasitic modes are responsible for the breakup of MRI 
modes into 3D turb ulence, although the exa ct role they play is 
still controversial dLesur & Longarettlll201 lb . Evidently in our 
ID model, parasitic instabilities are inhibited, allowing the pri- 
mary MRI mode to grow virtually for ever. To allow for the 
growth of parasitic modes, we have reproduced the same setup 
in 3D, adding 3D noise at moderate level ((v 2 ) ~ 10~ 2 ) to the 
ID perturbation described above. This systematically led to the 
production of an outflow before parasitic instabilities could do 
anything to the MRI eigenmodes. 

Although surprising at first sight, this result can be under- 
stood using the phenomenology of parasitic modes. First, it 
should be noted that the most uns table parasitic m odes are usu- 
ally Kelvin-Helmholtz modes (lLatter et al.l I2009D growing on 
MRI modes. The maximum growth rate of Kelvin-Helmholtz 
modes can be estimated by the local vertical shear rate. If we 
consider a primary MRI mode of amplitude 5v with a character- 
istic vertical size 61, then the maximum growth rate of the sec- 
ondary mode is y$ ~ Sv/6l. For the secondary mode to have 
an impact on the MRI mode, we require Vs > J w hich im- 
plies 6v > 5ly. Moreover, following Lat ter et aD (1201 Oh . we have 
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Fig. 1. Time evolution of the fiducial ID model IDRef. Left: magnetic field profile, right: Velocity and density profile. From top to 
bottom:? = 4.0; 7.0; 8.0; 60.0 
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Fig. 3. Streamlines (red dashed lines) and field lines (blue plain 
lines) of our steady solution obtained at t — 95. 



6b ~ Btfdv/QH where 6b is the magnetic perturbation of the 
MRI mode and B z q is the mean vertical field. Therefore, the par- 
asitic mode can destroy the MRI mode only if 



6b 61 y 
~Btf) > HQ. 



(10) 



In our system, the MRI mode characteristic length 61, growth 
rate y and the mean field amplitude are all of the order of 1 (in 
code units, see ^2. 11 1. This implies that parasitic modes will ap- 
pear only when 6b > 1. However, as it can be seen in fig. [1] 
the outflow starts when 6b < 1, i.e. before parasitic modes could 
act on the MRI mode. This result is due to the relatively large 
magnetisation used in these simulations (p < 1) compared to 
traditional MRI setups (p < 10~ 3 ) . This large magnetisation im- 
plies the production of a large scale MRI mode whose growth 
rate is of the order of the orbital time scale. This makes the out- 
flow more favourable compared to secondary instabilities. On 
the contrary, when the magnetisation is small, dominant MRI 
modes are found at smaller scale (large n, see fig.|2]l. In this case, 
parasitic modes are favoured and a turbulent flow is obtained. 

We should emphasize that the presence of an outflow does 
not mean that parasitic instabilities are totally absent from this 
picture. As we will show later (§@), solutions exhibiting an 
outflow are themself subject to parasitic instabilities. However, 
these instabilities have nothing in common with the traditional 
parasitic instabilities of MRI eigenmodes. 

We have seen above that the evolution of a large-scale MRI 
mode in a strongly magnetised shearing box leads naturally to 
the production of a magnetically-driven outflow. This steady 
outflow is essentially one-dimensional and can be described 
by v(z),p(z), B(z). In the following we will concentrate on the 
structure of this outflow: phenomenology, critical points, bound- 
ary conditions and conserved quantities. 

3.2. Outflow phenomenology 

As we have shown above, the outflow is primarily produced 
by the magnetic pressure gradient. The magnetic pressure be- 
ing maximum at z ~ 1.5, it pushes up the outflow at z > 2.0 
but it also compresses the bulk of the disc. An alternative view 
of this effect can be obtained looking at currents. The outflow 
is in this case due to horizontal currents which are reversed at 
z ~ 1.5. We typically have J x > and J y > in the bulk of the 
disc whereas J x < and J y < in the atmosphere z > 2. It is 
important to note that the outflow acceleration can occur only if 
these currents are non-zero and change sign. This remark justi- 
fies the absence of any outflow with "zero gradient" boundary 
conditions (see 93.41 1. 



It i s of interest to put t h ese o utflow solutions in the context 
of the Blandf ord & Pavnel ([1982) disc wind paradigm. In this 
model, the outflow is initiated by a magnetocentrifugal effect: 
the poloidal magnetic field lines are considered as rigid wires 
anchored in the bulk of the disc and fluid particles are allowed to 
drift along these wires. If the field lines are sufficiently inclined 
(typically more than 30° with respect to the vertical axis), then 
the particles are azimuthally accelerated by the anchored field 
lines. This leads to a centrifugal effect which ejects fluid parti- 
cles along field lines. In this picture, the ejection is driven by 
an exchange of angular momentum: angular momentum is taken 
from the disc by the field line and it is then released in the ejected 
material. 

In order to compare this mechanism to outflow solutions 
driven by MRI modes, we show in Fig. [3] the poloidal stream- 
lines and field lines of such a solution. We first note that poloidal 
streamlines and fieldlines are not aligned. This property is al- 
lowed by the shearing box boundary conditions, but it physi- 
cally means that magnetic flux is "accreted" toward the centre 
of the discfl In a more realistic model including curvature, such 
a state could not be sustained for a long time because magnetic 
flux would get accumulated at the disc centre, thereby modify- 
ing the disc properties (especially its rotational profile). This is 
the principle motivation for the presence of a strong "magnetic 
diffusivity" in disc wind models, either assuming the presence 
of small scal e turbulence ([Fe rreira & PelletierH l993a). ambipo- 
lar diffusion (IWardle & K6nigl ll993l) or Hall and Ohm diffusion 
dKonigl et al.ll2010l) . 

Despite this differenc e, we recover most of the p henomeno- 
logical properties of the Blandford & Pavne (1982) paradigm: 
field lines are inclined and drive an outflow which is inclined 
towards the same direction. This indicates that angular momen- 
tum is exchanged between the field and the flow. As we will 
see below ( jj3.5.31 l. angular momentum is effectively taken away 
from the disc by the field and then released to the ejected mate- 
rial. This effect is actually inevitable since the current configura- 
tion described previously inevitably leads to a positive magnetic 
torque oc J Z B X - J X B Z in the outflow (and a negative one in the 
disc midplane). Because angular momentum is taken from the 
disc by the field, a strong radial flow is produced which explains 
the streamlines' orientation forz < 1. Finally, we find an incli- 
nation angle of ~ 40° for the poloidal magnetic field lines and 
~ 25° for the poloidal velocity field. This last val ue is very close 
to the critical value of Blandfo rd & Pavnel (Il982l) . 

3.3. Critical points 

In principle, it is possible to look systematically for a steady ID 
solution of equations (Q} — (0. This is done writing the equations 
of motion in the form M • X = Y, where 



M = 



v z p ^ 

pv z -B z 

pv z -B z 

c 2 pv z B x B y 

-B, B x v z 

-B z B y v z 



(11) 



2 In a shearing box, the flux coming in the box through the x = +L v /2 
boundary condition is equal to the flux leaving the box through the 
x = -L x /2 boundary condition, making the magnetic field configura- 
tion overall stationary. Such a solution is however very specific to the 
shearing box and does not represent a steady situation in an accretion 
disc 
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Fig. 4. Vertical velocity and MHD wave speeds for the fiducial 
outflow solution at t — 60. 



X = d 7 
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2Qpv v 
-(2 - q)Qpv x 
-P&z 


-qO£ x 



(12) 



where M and Y are a matrix and a vector which do not contain 
any spatial derivative. In order to solve this nonlinear system, 
one then inverts the above system of equation to get a set of 
ordinary differential equations X = M _I ■ Y. However, when a 
critical point is reached, M is singular and the system cannot be 
inverted anymore. In this case, the physical system needs to sat- 
isfy an extra condition in order to get through the critical point. 
In the shearing box, we find three types of singular points: 



- two slow magnetosonic points 



z 2 



V 2 A + c 2 - A / ( y2 +c 2)2-4 c 2y2 ; 



- two Alfven points 



v. 



vi 



z - 'Az> 

and two fast magnetosonic points 
, 1 



Vj + c 2 + ^(V 2 + c 2 ) 2 - 4c 2 V 2 z 



(13) 



(14) 



(15) 



where Va = ~\jB 2 /p and Va z = B J -\/p. In the following, thanks 
to symmetries, we will only consider solutions with v. > so 
that only one critical point of each kind will be present. 

We present the MHD wave speeds and flow speeds in Fig. [4] 
We find that the slow point is located around z = 0.52 and the 
Alfven point is found at z = 2.47. The flow does not cross the fast 
point, however we find that the fast speed and the flow vertical 
speed tend to converge more rapidly close to the upper bound- 
ary. Note that a very similar behaviour was o bserved in global 
self-s imilar solution close to the fast surface (ICasse & Ferreira! 
2000, Fig. 1). This result indicates that the flow is still causally 
connected to the disc and therefore the boundary condition we 
impose at the top of the box still has an impact on the flow struc- 
ture itself. This point will be discussed in the next section. 



Run 


z B 




is 


z A 


Vb(Za)\ ) 


0v{za)( ) 


1 VnA 


A 


W. 14 / 


U.JO 3 


i.UJ 


JU.O 


o.u 


1Dz6 


6 


0.123 


0.530 


2.47 


58.6 


16.4 


1Dz8 


8 


0.105 


0.509 


2.94 


58.8 


22.1 


1Dz12 


12 


0.086 


0.483 


3.83 


58.7 


28.0 


1Dz16 


16 


0.072 


0.475 


4.77 


57.9 


31.4 


!Dz20 


20 


0.063 


0.470 


5.86 


56.5 


32.8 



Table 1. Evolution of the mass loss rate pv z , slow point zs and 
Alfven point za as a function of the altitude of the boundary 
condition zb- We also show the inclination angle with respect to 
z of the poloidal field lines and stream lines at the Alfven point 

(0b(za), 0v(za)) 



3.4. Influence of the upper boundary condition 

We have seen above that the outflow is still physically connected 
to the disc since it is not super-fast. Moreover, it looks as if the 
fast magnetosonic point is located close to the imposed upper 
boundary condition. One might wonder what is the exact role 
played by this boundary condition. In order to investigate this 
issue, we have performed two kinds of test, either modifying the 
altitude zb of the vertical boundary condition or modifying the 
nature of the upper boundary condition we apply. 

First, varying the altitude of the vertical boundary condition 
Zb led to the results presented in Tab.Q] One finds that modifying 
the altitude of the boundary condition strongly modifies the out- 
flow properties. In particular, we find that increasing the altitude 
of the boundary condition leads to a decrease in the flow ejec- 
tion rate. This evolution is accompanied by a slow point moving 
closer to the midplane and an Alfven point moving higher up 
in the atmosphere. This result clearly demonstrates that because 
the flow has not crossed the fast magnetosonic point, the solu- 
tion we obtain is still constrained by the boundary condition we 
impose at zb- In all the solutions described in Tab.Q] we have ob- 
served a "convergence" of the fast magnetosonic speed and the 
flow speed when approaching zb, similarly to what we find in 
Fig.0 This surprising result tends to indicate that our boundary 
conditions somehow force the fast point to be close to the upper 
boundary. 

The fact that the slow point moves closer to the disc mid- 
plane when the mass loss rate decreases might look dubious 
to readers familiar with the phenomenology usually associ- 
ated with the slow point. In particular, it is often said that the 
slow point "sets" the wind escape speed through the relation 
pVj. = Po(Zs)c wh ere po is taken as the hydrostatic density pro- 
file dSpruitlll996b . This argument predicts that smaller mass loss 
rate should be associated to slow points located higher up in the 
atmosphere, which is exactly the opposite of what we find. This 
simple argument is however not valid in the present case since 
the density profile deviates significantly from the hydrostatic 
equilibrium. In particular, configurations 1Dz16 and 1Dz20 ex- 
hibit strongly compressed discs due to a large magnetic pressure 
gradient in the atmosphere which most probably prevents signif- 
icant ejection from happening. 

We have also tried to modify the nature of the upper bound- 
ary conditions. First, instead of imposing B x (zb) = 0, we have 
imposed a fixed angle t o the poloidal fi eld, i.e. B x {zb) — tan(0)B z 
with 6 = 30° ; 45° as in lOgilvid d2012h . Surprisingly this did not 
modify significantly the outflow solution we obtained: the field 
values are modified by less than 5%. This can be explained by 
the fact that the outflow is super-Alfvenic when it reaches the 
top boundary. As we will see below ( ^3. 6i , sub-Alfvenic out- 
flows are effectively very sensitive to the field configuration at 
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the boundary, but super- Alfvenic outflows are not. We conclude 
from this that the inclination angle of the poloidal field line is 
set by the Alfven point crossing condition. This result is cor- 
roborated by the constant inclination angle at the Alfven surface 
found when one changes the box vertical size (Tab.Q]) 

We have finally tried to impose a zero gradient condition 
on B x and B y (classically called "outflow" boundary condition). 
This results in the suppression of the outflow solution. We ob- 
serve instead a constant increase of the magnetic pressure in 
the atmosphere which results in a strong compression of the 
disc material in the midplane until the disc occupies one nu- 
merical erid_c^lL_Thisj2sult is similar to the low B simula- 
tions of lHawlev et alJ (Q~995) with mean vertical flux. This was 
to be expected since the outflow is driven by horizontal currents. 
Imposing a zero gradient condition means that no horizontal cur- 
rent is allowed at the boundary, blocking any potential outflow 
by cancelling the vertical component of th e Lorentz force and 
the cha nge of sign of the magnetic torque dFerreira & Pelletierl 
Il993bh . 

3.5. Conserved quantities 

In the following, we assume all the quantities v, B, p 
depend only on z, as found in the steady ejection above. 
With this hypothesis, we reconstruc t the conserved quantities 
used in global disc wind models dBlandford & Pavnel [1982; 
iPelletier & Pudritzlll992l:ICasse & Ferre ira 2000). 



3.5.1. Frozen field condition 

We first note that under the above conditions, B- and pv z are con- 
stant in the box thanks to flux and mass conservation. We there- 
fore introduce a proportionality constant between these quanti- 
ties: 



3.5.2. Magnetic surface rotation 

A relation, known as the magnetic surface rotation speed can be 
deduced from the y component of the induction equation. Since 
By does not depend on x, we get after some algebra 



B P V 



—B v — v y + qQx 

P ' 



= 



(19) 



Note that the above relation is o nly valid along a po l oidal mag- 
netic field line. In contrast to Pelle tier & Pudrita (1 19921) . no 
equivalent relation is found along poloidal streamlines. This al- 
lows us to define a new conserved quantity 



Vy — qQ,X By 



(20) 



which is conserved along a magnetic field line. It should be 
pointed out at this stage that & x = u y B z - u z B y = B z v* where 
£> x is the x component of the electromotive force. Therefore, in 
the frame translating at v*, the x component of the electromotive 
force is zero, which justifies the naming of this invariant. 



3.5.3. Angular momentum conservation 

The angular momentum conservation equivalent is derived from 
the y component of the equation of motion (0. It should first be 
noted that this equation can be written: 



pv p -VZ-Bp- VB y = 



(21) 



where X. = v y + (2 - q)Q.x is the local equivalent of the spe- 
cific angular momentum. Since B x does not depend on x, we can 
rewrite the above equation using the flux freezing condition: 



pv p 



B 



(22) 



pv z = aB z 



(16) 



Thanks to the x induction equation, we also have v x B z - v z B x = 
S v = const, which can be simply written as: 



V r = Vr - V, 



El 

B- 



(17) 



where v* is a constant. Physically, v* is the advection speed of 
the poloidal magnetic field. In global models, this quantity is 
implicitly set to in order to avoid flux accumul ation at the disc 
inner edge and potential singularities at R = dChan drasekhar 
119561: lMeste]||1961h IPelletier & Pudritzl |1992|) . This is not re- 
quired here as we are using a shearing-box model. However, one 
should keep in mind that our model implicitly allows radial ad- 
vection of magnetic field lines. The existence of this conserved 
quantity allows us to define a relation between the poloidal field 
and the poloidal velocity, namely: 



a 

v p = -Bp 

P 



v,.e 3 



(18) 



where v p and B p are the poloidal (x, z) components of the ve- 
locity and magnetic fields. This equation shows a major differ- 
ence between the global approach and the local approach. In the 
global approach, v* = (no advection of magnetic field lines), 
which implies that poloidal streamlines and magnetic field lines 
are aligned. This is not the case in our solutions, for which 
v* = v x (z = 0) < 0. 



which indicates that / = £ - is conserved along streamlines. 

The angular momentum along one streamline is shown in 
Fig. [5] This allows us to check that the angular momentum is 
effectively conserved in our simulations. Moreover, we find that 
most of the angular momentum is initially extracted from the 
disc by the toroidal field B y . Higher above the disc (z ~ 2H), the 
angular momentum is exchanged between the toroidal field and 
the accelerated gas. This demonstrates that the magnetocentrifu- 
gal acceleration effect, present in the lBlandford & Pavnel (1982) 
phenomenology, is also found in our outflow solutions. 

3.5.4. Bernoulli invariant 

The Bernoulli invariant is derived from (0 after a scalar multi- 
plication by it. One gets: 



pu p 



—+H+rit 

2 V 



J X B -u 



(23) 



where H is the flow enthalpy. The right-hand side of this equa- 
tion corresponds to the work done by magnetic forces. This term 
can be calculated as a function of the previously defined invari- 
ants which, after some algebra leads to: 



pup ■ V 



+ H + ip- 



ByVl 



-qQ.pB y v* x . 



(24) 



Because v p and B p are not strictly speaking collinear, 
Bernoulli's equation is left with a non conservative term. 
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Fig. 5. Angular momentum along one streamline computed ac- 
cording to (1221 . The angular momentum is initially stored in the 
torooidal field before being transferred to the flow. 
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Fig. 6. Bernoulli invariant computed according to ( l24b . Ek'- 
Kinetic energy, Ej: thermal energy (enthalpy), Ef. Potential en- 
ergy, Eb'. Conservative magnetic energy, Eb,ic- non-conservative 
magnetic energy. 



However, as we will see later, this term is important only in the 
bulk of the disc, so that the flux defined above will be approxi- 
mately invariant along a streamline. Comparing this invariant to 
the one defined in global geometry, we note the presence of a 
new term in our case, B x v* x /a. As shown before, this term de- 
scribes the energy associated with the field lines being advected. 
We will see later the role it plays in the ejection process. 

We show the different terms involved in the Bernoulli invari- 
ant d24l) in Fig. [6] One should note that the magnetic contribution 
is separated in two parts: (i) the conservative part {B y v*-B x v x )/a 

and (ii) the non-conservative part J dl qQ.pB y v* x where the inte- 
gral is computed along a streamline. 

We first find that the invariant is approximately conserved for 
z > 0.2. Initially (up to z ~ 1.5) the energy is stored in the con- 
servative component of the magnetic field. The non-conservative 
part is constantly decreasing, demonstrating that this term is 
helping the outflow. Higher in the outflow, the magnetic energy 
is converted into kinetic and potential energy. Finally, we find 
that the thermal energy plays essentially no role in the ejection 
energetics. 

The fact that potential energy increases along a stream- 
line might look at first surprising since in a typical 

Blandford & Pavnel (Il982l) situation, one would expect the po- 
tential energy (gravitational+centrifugal) to decrease along a 
field line. This is not the case here because the inclination an- 
gle of the outflowing streamlines is slightly less than 30° (see 



3.6. Magnetisation dependency 

In the previous discussion, we kept a constant equivalent sur- 
face density X by artificially injecting mass in the disc midplane. 
This approximation, although partly motivated by the global 
disc structure, should be tested more extensively. To do so, we 
have performed simulations without mass injection. By defini- 
tion these simulations cannot achieve a steady state. They are 
however representative of an extreme case in which no mass is 
coming from the outer disc. 

It should first be pointed out that the outcome of these sim- 
ulations depends strongly on the nature of the upper boundary 
condition. This is because as the disc mass is lost, the Alfven 



Slow Point 

Alfven point 



'','iitS^ 



60 

t 



Fig. 7. Evolution of the critical points location as a function of 
time in a case without mass injection. The flow becomes sub- 
Alfvenic at t ^ 20 (see text). 



point moves higher up in the atmosphere (see Fig [7]). At some 
point, the Alfven point crosses the upper boundary and the ejec- 
tion becomes sub-Alfvenic. When this happens in this simpli- 
fied setting (no radial dependence), the poloidal field inclina- 
tion is not set anymore by the Alfven point crossing condition 
(Ogilvie 2012, see also 93.41 . Instead, it should be set manually 
at the upper boundary. If we set B x (zb) = 0, the ejection stops 
as soon as the Alfven point exits the simulation domain. This 
is expected since no ejection should be happening with strictly 
vertical poloidal field lines. On the contrary, if we impose a fixed 
inclination ang le of 45° at the upper boundary condition, as in 
OgilvT^ (1201 2l) the outflow is kept. In the following we will con- 
sider the latter case. 

Because the disc is losing mass, its magnetisation p is in- 
creasing as a function of time (see Fig.[9]l. We note however, that 
the system appears to be approaching a steady state with p ^ 2.3. 
Looking at the correlation between the mass loss rate due to the 
wind Mw = (pv z )z=z B anc ' me magnetisation p (Fig |7J, we find 
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Fig. 8. Flow magnetisation fi as a function of time in a case with- 
out mass injection. 




Fig. 9. Mass loss rate due to the wind as a function of the mag- 
netisation fi. Deduced from a case without mass injection 



that the mass loss rate decreases steeply for /i > 1, which ex- 
plains the quasi-stead y state we obser ve. We note that this last re- 
sult is very simila]0to Qgilvie> (l2012l) (Fig. 4). This demonstrates 
a plausible causal transition between MRI-drive n outflow solu- 
tions described here and the solutions found bv lQgilvid(l20T2h 
in the MRI stable regime. 

3. 7. Comparison with a global outflow solution 

Because the shearing box model is subject to several restric- 
tive hypothesis, it is useful to compare the MRI driven out- 
flows found here with global steady self-similar solutions of disc 
winds. The idea here is not to have a quantitive agreement be- 
tween these solutions since they are computed in very different 
different configurations. Instead, we look for qualitative similar- 
ities which could indicate that the physical processes at work 

3 Note that our definition of M w differs from lOgilvid d2012h by a 
factor V2tt. Moreo ver, our M w is normalized by the initial density in the 
midplane, whereas Ogilvie (2012) normalized M w by the instantaneous 
surface density, which in our case would decrease in time. 



are similar. T o this end, we have used the nume rical procedure 
describ ed by iFerreira & Pelletierl (Il995h and ICasse & Ferreiral 
d2000h to compute an isothermal "cold" solution with the fol- 
lowing parameter^ e — 0.1 ; a m — 1.0 ; Pm — ; fi — 0.7268 ; 
£, = 9.92x 10~ 3 . This last parameter is connected to the radial de- 
pendence of the accretion rate through the relation M a (r) cc . 
Although this solution depends explicitly on r and z through a 
self-similar ansatz, the radial dependence is much weaker than 
the vertical one up to a few scale heights. Therefore, we only 
show the z dependence of the resulting solution in Fig. \W\ We 
have zoomed on the region from z = . . . 6H to allow a com- 
parison with shearing box solutions (Fig.[TJ, although the actual 
global solution extends up to 25H (fast magnetosonic point). 

The direct comparison between this global solution and the 
shearing box solutions shows several differences: the velocities, 
and in particular v z , are larger in the global solution. On the other 
hand, the horizontal magnetic field strength is weaker compared 
to the shearing box solution. However, we recover most of the 
physical properties found in global solutions in the shearing box 
solution: a strong magnetic shear close to the midplane accom- 
panied by a magnetically compressed disc (compared to a purely 
hydrostatic equilibrium) and a strong accretion flow in the disc 
midplane (v x 2c). 

It should be noted that the global solution presented here has 
a much weaker mass loss rate (pv z ~ 5 x 10~ 3 ) compared to 
the shearing box solution (pv z 0.1). This can be partly due to 
the stronger mean vertical field required by the global solution 
(fi ~ 0.24), although this difference is not enough to explain 
entirely this discrepancy (Fig.|9]l. The fact that the mass loading 
is much smaller in the global solution explains the faster escape 
velocity found in this solution. 

4. Three dimensional solution and stability 

4. 1 . Outflow evolution in 3D 

In order to investigate the 3D evolution of the ID outflow desr- 
cribed in §3. II we have perform the outflow simulations in 
3D. The initial condition consists of the ID initial perturba- 
tion described in ||3]plus a random 3D noise added to v x with 
max(|v x |) = 0.02. We show in Fig.[TT]the temporal evolution of 
simulation 3DRef for /j = 8x 10~ 2 . These figures represent the 
evolution of the density, poloidal magnetic field inclination and 
vertical velocity averaged in the (x, y) plane in a (z, t) diagram. 
We first observe the presence of a strong outburst (t ~ 8) fol- 
lowed by a rather steady state during which the flow does not 
evolve rapidly. This state corresponds to the ID solution de- 
scribed in ^3] and is essentially a ID outflow solution. However, 
after some time, this ID outflow goes unstable (f ~ 40) and rapid 
variations in all quantities are observed. Surprisingly, the global 
structure of the ID outflow is conserved: on time and horizontal 
averages, the typical v z vertical and poloidal inclination profiles 
are consistent with the ID steady solution (Fig. [TJi. We therefore 
produc e a turbulent outflow driven by the Blandfor d & Pavnel 
(ll982l) magnetocentrifugal effect. 

We show in Fig. [12] snapshots of the same simulation taken 
at t = 20, t = 40 and t = 270. At t = 20, we confirm the ID 
nature of the solution: no dependence can be seen in x or y. This 
also shows that the outflow is produced before parasitic modes 
get significantly excited. At t = 40 we see the development of a 



4 The precise defini tion of these parameters can be found in 
ICasse & Ferreiral i2000l) . Note that the definition of the fi parameter 
used in these solutions differs slightly from the /i used in this paper. 



9 



Geoffroy Lesur, Jonathan Ferreira and Gordon I. Ogilvie: The magnetorotational instability as a jet launching mechanism 




Fig. 10. Global solution computed following Ferreira & Pelletierl d!995h in the z = . . . 6H region (see text). Colours and units 
identical to Fig.Q] 



"secondary" instability of the outflow solution which produces 
a ^-dependent structure. The physical processes responsible for 
this instability will be discussed in §14.21 When the "secondary" 
modes reach a significant amplitude, the structures break down 
into non axisymmetric turbulent motions in which the main out- 
flow properties are maintained (field line inclination, outflow 
speed). A typical example of such a state is shown at t = 270. 

4.2. Outflow solution stability 

We have shown above that the outflow solution is unstable 
to 3D perturbation, resulting in a "turbulent outflow" config- 
uration. Several authors have discussed the possibility of hav- 
ing such an instability (Lubow et all 11994 : Lubow & Spruit 
119951; ICao & Spruitl 120021) although the applicability of these 
stability analyses to a l l outflow so lutions is still uncertain 
(iKonigl & Wardldll99rtlKo^l2004h . 

In order to analyse the instability observed in our outflow so- 
lutions, we have performed a simulation (3DLin) starting from 
the ID solution corresponding to the final state of run IDRef to 
which a small-amplitude (10 -3 ) 3D white noise was added. Since 
the growth phase of the instability implies only jc-dependent 
modes, we use a Fourier decomposition in the x direction to 
characterise growing modes: 



p = po(z, t) + 5p(x, z, t) 
6p(x,z,t) = !R [ ^ pk{z, t) expplx) 



(25) 
(26) 



and similarly for v and B. In this expansion, we have assumed 
the instability was growing on the top of the ID solution po(z) 
given by the final state of run IDRef. We first present the tem- 
poral evolution of the fluctuation ^j{5p 2 ) in Fig.Q~3] We find that 
perturbations grow exponentially with a growth rate y = 0.33Q 
up to t ^ 40 where a saturation is reached. To further inves- 
tigate the behaviour of this instability, we present the temporal 
evolution of \pk(z = 0.5)| as a function of k in Fig. [14] As it 
can readily be seen, the growth is dominated by modes having 
n = kL x /2ji = 4 which is consistent with the 4 "spots" ob- 
served in Fig. [T2] The measured growth rate of this mode is 
y = 0.33f2 which explains its fast appearance in 3D simula- 
tions once an outflow has formed. However, other neighbouring 
modes are growing as well (n = 3; 5; 6) although not as fast as 
the n — 4 mode. Finally once the n — 4 mode reaches large am- 
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Fig. 12. Snapshots of the p = 8 x 10~ 2 solution. Logarithm of 
the density is represented in coloured volume rendering whereas 
field lines are represented as tubes. From top to bottom: t -20; 
t = 40 ; t = 270. x axis is horizontal and z axis is vertical. 



plitudes (|pa | ~ 1), we note the sudden growth of the n — 8 and 
n — 9 modes which are the result of nonlinear interactions of the 
fastest growing modes n — 3; 4; 5; 6. 

In order to analyse the physics underlying the n — 4 mode, 
we present in Fig. [15] the density fluctuations corresponding 
to that eigenmode. We first note that the density fluctuation is 
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Fig. 11. Spacetime diagram of horizontally averaged quantities in afi — 8x 10 2 simulation. Top: log 1() (p), middle: poloidal magnetic 
field line inclination with respect to the vertical axis(in °), bottom: vertical velocity. 



highly localised in z around z ~ 0.5. In addition, this eigenmode 
has a tail whose inclination and shape closely follow that of the 
poloidal magnetic field (Fig. [3). The localisation of the eigen- 
mode is rather surprising and requires some explanation. We first 
note that this localisation is much higher than the top of the mass 
injection region {H m ^ = 0.1). However, comparing the relative 
vorticity component a> v = d : v x of the ID outflow (Fig. [T6l l to 
the vertical profile of the eigenmode (Fig [TTb . we find that the 
density perturbation is localised close to a maximum of u) x . This 
tends to suggest that this instability is somehow linked to the ver- 
tical (L> y profile of the outflow solution, and therefore to a kind of 
Kelvin-Helmholtz instability. This potential link is also consis- 
tent with the growth rate y < max(w v ) and with the shape of the 
eigenmode and the vorticity profile around z = 0.5. 

We would like to stress that these remarks are not a proof that 
this outflow instability is of the Kelvin-Helmholtz type. Among 
the effects we did not take into account in that analysis are the 
magnetic field, compressibility and the presence of a large v z up 
in the atmosphere. However, if we assume that the source of the 
instability lies around z ~ 0.5 as suggested by the eigenmodes, 
then time scales, length scales and phenomenology match that 
of a Kelvin-Helmholtz instability. To ascertain these claims, a 



proper linear study taking into account all of the outflow proper- 
ties is required, which is well beyond the scope of this paper. 



5. Discussion 

In this paper, we have shown that large scale MRI modes which 
are unstable when the disc magnetisation is moderately sub- 
thermal spontaneously produce super- Alfvenic outflows. The 
physical mechanism behind this outflow is a Blandford & Payne- 
like process where angular momentum is transferred to bending 
field line and then released to accelerated material. We demon- 
strated that this outflow is qualitatively similar to the outflow 
solutions found both in local boxes and in global self-similar ge- 
ometry, making a clear connection between the MRI and the for- 
mation of disc winds. We have also shown that MRI outflows are 
unstable in 3D which could be a potential source of variability 
for disc winds and jets. These 3D instabiliti es could also be the 
origin of the turbulent resistivity a m used in lFerreira & Pelletierl 
( 1995) and subsequent self-similar models. 

However, the picture provided here is still incomplete. We 
first note that the simulations we have performed here only 
produce super-Alfvenic flows, whereas real escaping outflows 
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Fig. 13. RMS amplitude of the fluctuations 6p in run 3DLin. We Fig. 16. y component of the vorticity in the stationary outflow 
observe a linear growth phase with y = 0.33£1 solution. 
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Fig. 14. Amplitude of the Fourier modes \pk(z = 0.5)| during the 
linear phase of the outflow instability. Mode decomposition are 
plotted every A? = 2 from t = to t = 40. The instability is 
dominated by the n — 4 mode (see text). 




Fig. 15. Density fluctuations corresponding to the n — 4 eigen- 
mode. 



should be faster than the fast magnetosonic speed. This problem 
seems to be linked to the shearing box geome try, as several other 
autho r s have noticed this in shearing boxes dSuzuki & Inutsuka 
2009, Froman g" et alJ 1201 2l X. Bai private communication). 
Indeed, a shearing box does not allow for the opening of mag- 
netic field lines which is expected in realistic global geometry. 
We believe that such an opening is required to get super-fast 



Fig. 17. Eigenmode n — 4 vertical density profile 



flows, which are therefore beyond the scope of our shearing 
box model. The fact that the flow is only super-Alfvenic indi- 
cates that the upper boundary condition still plays a role in de- 
termining the outflow structure, and indeed it does, at least for 
the mass-loss rate ( 133.41 ). 

We also emphasize that the shearing box model pos- 
sesses a horizontal symmetry by the transformation (v x , v y ) — > 
(-v x ,-v y ), (B x ,B y ) -> (-B x , -By) and (x,y) -> (-x,-y). This 
symmetry indicates that locally, there is no mathematical differ- 
ence between a magneto-centrifugally accelerated wind where 
angular momentum is transferred from the disc to the jet and a 
magnetically decelerated accretion column (formally a magneto- 
centripetal wind) where the angular momentum is transferred 
from material falling radially inward to the disc. This symme- 
try is obviously broken once curvature terms are taken into ac- 
count. The presence of this symmetry implies that shearing box 
simulations can spontaneously switch from one situation to the 
other, which is unexpected in realistic situations. Such sudden 
changes in the magnetic configuration w ere indeed observed i n 
rare occasions in our 3D runs but also by Froma ng et alJ d2012l) . 

Finally, we should point out that the presence of a non-zero 
toroidal electromotive force implies that magnetic field lines are 
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accreted. As mentioned earlier, this situation is rather unrealis- 
tic in global geometry, although it is allowed in shearing boxes. 
More realistic configurations, probably including a sort of resis- 
tivity (either effective or molecular) are required to compensate 
for this effect. 

All of these points indicate that shearing boxes are not very 
well suited to study globally the outflows produced by MRI tur- 
bulent accretion discs. In particular, little can learned regarding 
outflow mass loss rate and velocities. We note however that our 
solutions can be qualitatively compared to global solution and 
several properties are recovered by the local model. Moreover, 
3D instabilities can be studied much more easily in boxes than 
in global geometry where computational costs increase very 
rapidly. 

Overall, our results tend to indicate a paradigm shift: up to 
now, the MRI driving "viscous" discs and disc winds at the ori- 
gin of jets have often been considered as separate processes. 
Here we show that these two processes are actually intrinsically 
connected: outflows are a logical consequence of the MRI in 
strongly magnetised discs. Obviously, the next question is to un- 
derstand what is driving the disc magnetisation. This could be 
due to local dynamos, large scale field redistribution through ad- 
vection and diffusion, etc. To identify these processes, shearing 
boxes are clearly insufficient and global models, including both 
large-scale magnetic fields and turbulence, will be required. 
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Run 


Resolution 


Box size 


Mass injection 


Smoothed potential 


Boundary condition 


Outflow 


IDRef 


1 x 1 x 256 


1.0 x 1.0 x 6.0 


Yes 


No 


BJzr) = 


Yes 


1Dz4 


1 X 1 x 170 


1.0 x 1.0 x 4.0 


Yes 


Yes 


BJzr) = 


Yes 


1Dz6 


1 X 1 x 256 


1.0 x 1.0 x 6.0 


Yes 


Yes 


BAzr) = 

XK^B / " 


Yes 


1Dz8 


1 x 1 x 340 


1.0 x 1.0x8.0 


Yes 


Yes 


BAzb) = 


Yes 


1Dz12 


1x1x512 


1.0 x 1.0 x 12.0 


Yes 


Yes 


B x {z B ) = 


Yes 


1Dz16 


1 x 1 x 684 


1.0 x 1.0 x 16.0 


Yes 


Yes 


BAzb) = 


Yes 


1Dz20 


1 x 1 x 856 


1.0 x 1.0x20.0 


Yes 


Yes 


BAzb) = 


Yes 


1DZG 


1 x 1 x 256 


1.0 x 1.0x6.0 


Yes 


No 


d-BAzs) = 


No 


1 DTwr 




1 n v 1 n v 6 n 


ICS 




J BAzb) = B z (zb) 
\ BAzb) = 


ICS 


IDNoMass 


1 x 1 x 256 


1.0 x 1.0x6.0 


No 


No 


BAzb) = 


No" 


IDNoMassInc 


1 x 1 x 256 


1.0 x 1.0x6.0 


No 


No 


f BAzb) = BAzb) 
\ B y (z B ) = 


Yes 


3DRef 


128 x 128x256 


16.0 x 16.0x6.0 


Yes 


No 


BAzb) = 


Yes 


3DLin" 


128 x 128x256 


16.0 x 16.0x6.0 


Yes 


No 


BAzb) = 


Yes 



Table 2. List of the simulations discussed in this paper. 

": The outflow disappears when the Alfven point gets out of the simulation box ( 33.61 ) 

b : The initial condition for this simulation corresponds to the final state of IDRef to which small amplitude 3D white noise was 
added ( $47% 
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